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Abstract. - We report Monte Carlo results for the two-dimensional hard disk system. Simula- 
tions were performed in the NVT ensemble with up to 65536 disks, using a new updating scheme. 
We analyze the bond orientational order parameter and correlation length in the isotropic phase 
and the scaling behaviour of the bond orientational order parameter in the transition region. The 
data are consistent with predictions of the Kosterlitz-Thouless-Halperin-Nelson- Young theory, 
while a first-order phase transition is unlikely and a one-stage continuous transition can be ruled 
out. 



The nature of the two-dimensional melting transition has been an unsolved problem since 
many years [[[J [|. The Kosterlitz-Thouless-Halperin-Nelson- Young (KTHNY) theory || pre- 
dicts two continuous transitions. The first transition occurs when the solid (quasi-long-range 
positional order, long-range orientational order) undergoes a dislocation unbinding transition 
to the hexatic phase (short-range positional order, quasi-long-range orientational order). The 
second transition is the disclination unbinding transition which transforms this hexatic phase 
into an isotropic phase (short-range positional and orientational order). There are several 
other theoretical approaches for the transition. One alternative scenario has been proposed by 
Chui 0). He presented a theory via spontaneous generation of grain boundaries, i.e. collective 
excitations of dislocations, and predicted a conventional first-order phase transition from the 
solid to the isotropic phase. In this case there exists a region where both phases coexist instead 
of a hexatic phase. Even for the simple hard disk system no consensus about the existence of 
a hexatic phase has been established. 

The melting transition of the hard disk system was first seen in a computer simulation by 
Alder and Wainwright ||. They used a system of 870 disks and molecular dynamics methods 
(constant volume V, energy E and number of particles N simulations) and found that this 
system undergoes a first-order phase transition. But the results of such small systems are 
affected by large finite-size effects. Recent simulations used Monte Carlo (MC) techniques 
either with constant volume (NVT ensemble) || fjj or constant pressure (NpT ensemble) 
[|[ Lee and Strandburg || used isobaric MC simulations and found a double-peaked 

Typeset using EURO-TEX 



2 



EUROPHYSICS LETTERS 



structure in the volume distribution. Lee-Kosterlitz scaling led them to conclude that the 
phase transition is of first order. However, the data are not in the scaling region, since their 
largest system contained only 400 particles. MC investigations of the bond orientational order 
parameter via finite-size scaling with the block analysis technique of 16384 particle systems 
were done by Weber, Marx and Binder [Q. They also favoured a first-order phase transition. 
In contrast to this, Fernandez, Alonso and Stankiewicz |J predicted a one-stage continuous 
melting transition, i.e. a scenario without a hexatic phase. Their conclusions were based 
on the examination of the bond orientational order parameter in very long runs of different 
systems up to 15876 particles and hard-crystalline wall boundary conditions. The analysis of 
Zollweg and Chester || for the pressure gave an upper limit for a first-order phase transition, 
but is compatible with all other scenarios. 

In this letter, we present results obtained through MC simulations in the NVT ensemble 
to answer the question of the kind of the phase transition. We consider systems of N = 32 2 , 
64 2 , 128 2 and 256 2 hard disks in a two-dimensional square box. We find that finite-size effects 
with these boundary conditions are not substantially larger than in a rectangular box with 
ratio y/3 : 2, furthermore no simulations in the solid phase were made. The disk diameter is 
set equal to one. For the simulations a new updating scheme was developed fl(i|| , in which the 
conventional Metropolis step of a single particle is replaced by a collective (non-local) step of 
a chain of particles. A cell structure was chosen such that one cell can only be occupied by 
a single disk. In all updatings, the random number generator proposed by Liischer [p"l| was 
applied. The simulations were performed on a Silicon Graphics workstation and a CRAY T3E. 
In the latter case, we used the different nodes of the parallel machine to generate independent 
data sets. Statistical errors have been calculated by binning. Additionally, we performed a 
jackknife analysis of the different data sets from the different nodes. Careful attention has 
been paid to the equilibration of all systems. For example, we performed 3 x 10 5 'sweeps' for 
N = 256 2 at p — 0.890 with the improved (chain) Metropolis updating scheme to warm up 
the system and 1.9 x 10 6 'sweeps' to measure the expectation values (for 6 independent data 
sets). The acceptance rate for this run was about 54%. Further details will be published later 



Simulations were performed in the isotropic phase and in the phase transition region. In 
the isotropic phase we measured the (global) bond orientational order parameter ^6 and the 
correlation length of the bond orientation ^6- The local value of ip6 for a particle i located at 
x = (x, y) is given by 



where the sum on j is over the iVj neighbours of this particle and 9ij is the angle between the 
particles i and j and an arbitrary but fixed reference axis. Neighbours are obtained in a usual 
way by the Voronoi construction. The (global) bond orientational order parameter is just the 
absolute value of the average over all particles: 



The bond orientational correlation length was extracted from the 'zero-momentum' correlation 
function of ip& (x) 
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by fitting the data with a single cosh, where Nk denotes the number of particles in a stripe 
between x + Ax/2 and x — Ax/2. This method allows a precise measurement of ipe apart from 
some systematical errors, which will be large — compared to our statistical errors — for small 
values of (for details sec [|l2]). Additionally, we calculated the radial bond orientational 
correlation function 

96(r) = {i> 6 *(a)Mr))/9(r) , (4) 

and extracted the correlation length from an Ansatz of the form g^r) ~ r -1 / 2 exp(— r/£ 6 ), 
where g(r) is the pair correlation function. In all simulations of the isotropic phase we used 
systems of at least TV = 64 2 particles and chose the box length to satisfy L > 7^6 . In 
these cases, within statistical errors, we found no finite-size effects on the bond orientational 
correlation length and on the susceptibility 

X6 = N(^ 6 2 ) . (5) 

Equation (|5|) differs from \e = N((ip e 2 ) — (i/Je) 2 ) by a factor 1 — 2/tt in the thermodynamic 
limit. Due to the new updating scheme, which reduces the autocorrelation time, we were able 
to perform simulations with large correlation lengths, namely close to the disclination binding 
transition point p\. 

The KTHNY scenario predicts an exponential singularity for the correlation length 

&(*)=a 4 expfV- 1 / 2 ) , (6) 

and the susceptibility X6 

Xe(t) = a x exp^i- 1 / 2 ) (7) 
if t = pi — p — > + . The critical exponent rj e defined by 

Xe ~ & 2 -" 6 , (8) 
is given by 776 = 1/4, while is a non-universal constant and 

b x = (2 - % )6 C . (9) 



Wc analyzed the critical behaviour of £g and \ e by performing least square fits according 
to eqs. (||) and (0). Typically, statistical errors of £ 6 an( l Xe are 111 the range 1% — 5%. 
Errors for the fitting parameter were computed by performing fits on data sets being Gaussian 
distributed around the expectation value. If we used all 12 different measurement points 
with 0.82 < p < 0.89 we got a x 2 P er degree of freedom (d.o.f.) of 0.75 for £g(t) and 0.65 
for X6 (t) , i- e. the data are in a very good agreement with an exponential singularity of the 
KTHNY type. The critical values of p were given by pi = 0.9017(7) and 0.9002(3), respectively. 
The results for the susceptibility are shown in fig. [j]. Data far away from the transition 
(in particular for the correlation length) are affected by systematical errors. Therefore, fits 
were also performed omitting some data at low densities. For example, for the eight points 
with 0.855 < p < 0.89 (for £ 6 > 3.0) we got x 2 /d-o.f.= 0.23, p K = 0.9006(8) for £ 6 (i) and 
X 2 /d.o.f.= 0.58, pi = 0.9000(4) for Xe(t)- The critical exponent rj 6 is calculated using eq. ([)[), 
yielding t]q — 0.451(21) and 0.349(44), respectively. A detailed analysis shows that the value 
of ?]6 decreases if t — > + . This can be seen, if we use eq. (ft) and plot hi(x6/£6 7 ^ 4 ) versus 
In(£e)- For the predicted value = 1/4 we should see a horizontal line. A different value of 
776 would correspond to a straight line with a non-zero slope. Indeed, there is a negative slope 
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Fig. 1. - Susceptibility as a function of the density. The curve shown is the best fit for a KTHNY 
behaviour. The critical value of p is visualized by a vertical line. 
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Fig. 2. - Test of the scaling relation xe ~ £.1^- The slope gives the deviation from 776 = 1/4. 



for small values of ^6 a s can be seen in fig. g. Nevertheless, in the limit ^6 — > 00 the data are 
compatible with rj 6 = 1/4. A fit with the last six data points gives % — 0.251(36). 

We now come to the simulations with p « p\. Finite-size scaling (FSS) implies X6 ~ L 2 ~ riB 
at p — pi for large enough systems. For p < p u corrections for finite correlation lengths of 
0(L/£g) have to be taken into account. For p- x < p < p m , t] 6 is a function of the density. As 
p approaches the melting density p m , i.e. at the end of the hexatic phase, rje —> + . Figure 
I shows ln(x 6 /i 7/4 ) versus ln(L) for various p. The slope was extracted from linear fits and 
gives the deviation from rje = 1/4. Using the FSS behaviour to locate p u the requirement 
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Fig. 3. - Finite-size scaling of the susceptibility for various densities. The slope gives the deviation 
from r?e = 1/4. The lines are guides to the eye. 

VeiPi) = 1/4 yields p\ — 0.899(1). This value is in agreement with that obtained from the 
singularities of £e(t) and X&if)- A slightly different value of rje would not alter this situation. 
Moreover, our estimates of pi agree with Weber et al. (?J who used the fourth-order cumulant 
intersection (p^ — 0.8985(5)). However, it differs from their value obtained using the singularity 
of X6 {fix — 0.913). The result p\ = 0.916(4) of Fernandez et al. || is not compatible with our 
value. 

Our MC data suggest that the the orientational order behaves as predicted by the KTHNY 
theory. A one-stage continuous transition || can be ruled out, since p m > 0.910 (obtained 
from i](p m ) = 0) is away from p\ in this work. Also our data are not compatible with a 
first-order phase transition with small correlation lengths, since we find no deviation from the 
predicted singularities of \6 and ^6 up to ^6 ~ 38. (Alternative approaches for the singularities 
result in large x 2 /d.o.f. For example, a conventional second-order behaviour with a power-law 
singularity of the form ln(£ 6 ) = a — v\n(t) yields x 2 /d.o.f.=4.1.) We also have examined FSS 
of the fourth-order cumulant U = 1 — (ipe ) /3(ipQ 2 ) 2 . In the hexatic phase, FSS implies scale 
invariance of U. The narrowness of such a scale invariant region was one argument in ref. 
against the existence of a hexatic phase. Unfortunately, statistical errors in our data are too 
large to answer this question. Although our data cannot rule out a first-order phase transition 
with very large orientational order correlation lengths, a KTHNY-like phase transition seems 
to be more likely. 
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